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Abstract 

The  paper  demonstrates  how  large-scale  stochastic  linear  programs  with  re¬ 
course  can  be  efficiently  solved  by  using  a  blending  of  classical  Benders  decomposi¬ 
tion  with  a  relatively  new  technique  called  importance  sampling.  Numerical  results 
of  large-scale  problems  in  the  area  of  expansion  planning  of  power  systems  and 
financial  planning  are  presented. 


1.  Introduction 

Solutions  obtained  from  deterministic  planning  models  are  usually  unsatisfac¬ 
tory  because  they  fail  to  hedge  against  unfavorable  events  which  may  occur  in  the 
future.  Stochastic  models  address  this  shortcoming,  but  in  the  past  have  seemed  to 
be  untractable  because,  even  for  a  relatively  small  number  of  parameters,  subject 
to  uncertainty  the  size  of  such  problems  can  get  very  large.  Stochastic  problems 
have  been  studied  extensively  in  the  literature  since  Dantzig  (1955),  for  example 
Birge  (1985),  Ermoliev  (1983),  Frauendorfer  (1988),  Higle  and  Sen  (1989),  Kail 
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fice  of  Naval  Research  Contract  N00014-89-J-1659;  the  National  Science  Foundation 
Grants  ECS-8906260,  DMS-8913089,  the  Electric  Power  Research  Institute  Con¬ 
tract  RP  8010-09,  CSA-4O05335,  and  the  Austrian  Science  Foundation,  “Fonds  zur 
Forderung  der  wissenschaft lichen  Forschung,”  Grant  J0323-Phy. 

Any  opinions,  findings,  and  conclusions  or  recommendations  expressed  in  this 
publication  are  those  of  the  authors  and  do  NOT  necessarily  reflect  the  views  of  the 
above  sponsors. 
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(1979),  Pereira  et  al.  (1989),  Rockafellar  and  Wets  (1989),  Ruszczynski  (1986), 
Wets  (1984)  and  others  contributed  in  this  area  (Ennoliev  and  Wets  (1988)).  The 
paper  addresses  a  common  class  of  stochastic  models,  namely  the  stochastic  linear 
programs  with  recourse.  We  discuss  how  decomposition  techniques  and  importance 
sampling  can  be  applied  to  solve  them  (Dantzig  and  Glynn  (1990),  Infanger  (1990)). 

2.  A  Class  of  Multi-Stage  Stochastic  Linear  Programs 

Large-scale  deterministic  mathematical  programs,  used  for  operations  and 
strategic  planning,  often  are  dynamic  linear  programs.  These  problems  have  a  stair¬ 
case  (multistage)  matrix  structure.  In  general,  the  size  of  these  stochastic  problems 
can  get  extremly  large  because  the  number  of  scenarios  grows  exponentially  with 
the  number  of  periods.  We  will,  however,  address  a  certain  restricted  class  whose 
number  of  scenarios  grows  linearly  with  the  number  of  stages:  The  problem  (whose 
constraints  are  stated  below)  breaks  down  into  two  parts:  a  deterministic  dynamic 
part  and  a  stochastic  part.  We  call  the  deterministic  part  the  master-problem.  It 
is  a  dynamic  linear  program  with  T  stages.  The  vectors  ct  and  6*,  and  the  matrices 
-Bt_i  and  At  are  assumed  to  be  known  with  certainty. 

T  T 

min  £<*  I,  +£*W) 

t=i  <=i 

—Bt-iXt-i  +  At  Xt  =  bt,  t  =  1, . . . ,  T,  Bo  =  0 

xt  +Dty'j't  =  d?',  t  =  1, . . .  ,T,  ut  €  ft* 

V?  >  0 

Each  stage  is  associated  with  a  stochastic  sub-problem.  Uncertainty  appears 
in  the  recourse-matrix  Ff*  and  in  the  right  hand  side  vector  where  u>(,  denotes 
an  outcome  of  the  stochastic  parameters  in  period  t,  with  denoting  the  set  of  all 
possible  outcomes  in  period  t.  The  sub-problems  in  each  stage  are  assumed  to  be 
stochastically  independent.  The  sub-problem  costs  ft  and  the  technology  matrix 
Dt  are  assumed  to  be  deterministic  parameters. 


2 


Facility  expansion  planning  is  an  example  of  this  type  of  formulation.  The 
master-problem  models  the  expansion  of  the  facilities  over  time.  Decision  variables 
are  the  capacity  built  and  the  capacity  available  at  time  t.  The  sub-problems  model 
the  operation  of  these  capacities  in  an  uncertain  environment.  Take  for  example 
the  case  of  expansion  planning  of  power  systems:  The  expansion  or  replacement  of 
capacities  of  generators  and  transmission  lines  are  determined  in  the  master  prob¬ 
lem.  The  capacities  at  each  period  t  are  made  available  to  the  system  for  operation. 
The  subproblems  model  the  power  system  operation,  the  optimal  sceduling  of  the 
available  capacities  to  meet  the  demand  for  electricity.  The  availabilities  of  gener¬ 
ators  and  transmission  lines  and  the  demands  are  uncertain  and  not  known  at  the 
time  when  the  expansion  decision  is  made. 

The  approach  is  primarily  “here  and  now”  (Dantzig  and  Madansky(1961)) 
and  justified  by  high  investment  costs  and  long  lead-times  for  capacity  expansion. 
However,  as  the  operations  subproblems  are  stochastically  independent  and  only 
expected  operation  cost  rather  than  the  state  of  the  system  after  period  t  affects 
the  expansion  plan  (as  failures  of  equipment  get  repaired,  and  uncertainty  in  the 
demands  are  interpreted  as  deviations  from  a  demand  path),  ’’here  and  now”  is 
equivalent  to  “wait  and  see”.  That  means  that  the  optimal  decision  in  period 
t  +  1  depends  only  on  the  capital  stock  on  hand  at  the  start  of  period  t  +  1  and  is 
independent  of  any  observed  outcomes  in  period  t,  i.e.  the  same  optimal  capacity 
expansion  decision  would  be  made  before  as  after  period  t  operations.  Thus  the 
facility  expansion  plan  can  be  laid  out  at  the  beginning  for  the  whole  planning 
horizon  based  on  the  expansion  costs  and  the  expected  operation  costs.  This  permits 
the  multi-stage  problem  to  be  treated  as  if  it  were  a  two-stage  problem.  The  first 
“stage”  concerns  the  single  decision  of  what  facility  expansion  will  be  in  all  future 
periods  without  knowledge  of  the  particular  outcomes  of  the  uncertainty  parameters 
in  future  periods.  The  second  “stage”  concerns  the  operations  problems,  where  the 
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recourse  decisions  made  depend  cm  the  realizations  of  the  stochastic  parameters. 
Note  that  for  T  =  1,  the  problem  is  exactly  a  two  stage  stochastic  linear  program 
with  recourse.  For  T  >  2  the  problem  is  a  two  “stage”  problem  with  the  second 
stage  consisting  of  T  independent  subproblems. 

3.  Multidimensional  Integration 

The  difficulty  of  solving  large-scale  stochastic  problems  arises  from  the  need  to 
compute  multiple  integrals  or  multiple  sums.  The  expected  value  of  the  second  stage 
costs  in  period  t  (we  suppress  the  index  t  for  this  discussion),  e.g.  z  =  E(fyu)  = 
E(C)  is  an  expectation  of  functions  C(vu)tw  6  ft,  where  C(t/W)  is  obtained  by  solv¬ 
ing  a  linear  problem.  V  (in  general)  is  a  /i-dimensional  random  vector  parameter, 
e.g.  V  =  (Vi, . . . ,  Vfc),  with  outcomes  vu  =  (t>i, . . .  ,Vh)u-  For  example  Vi  represents 
the  percent  of  generators  of  type  i  down  for  repair  or  transmission  lines  not  operating 
and  v?  the  observed  random  percent  outcome.  The  vector  «"  is  also  denoted  by  v, 
and  p(v“)  alias  p(v)  denote  the  corresponding  probability,  ft  is  the  set  of  all  possible 
random  events  and  is  constructed  by  crossing  the  sets  of  outcomes  ftt-,  *  =  1, . . . ,  h  as 
ft  =  fti  xft2  X  •  •  -xftfc.  With  P  being  the  probability  measure  under  the  assumption 
of  independence  the  integral  E  C{V)  =  /  C(vw)P(doj)  takes  the  form  of  a  multiple 
integral  E  C(V)  =  /♦•♦/  C(v)p(v)dv i . . .  dvh,  or,  in  case  of  discrete  distributions, 
the  form  of  a  multiple  sum  E  C(V)  =  C'Ct^Pi^i) . .  .p*(vfc). 

In  the  following  discussion  we  concentrate  on  discrete  distributions.  This  is  not 
a  restriction  as  the  approach  can  be  easily  adapted  for  continuous  distributions.  In 
practical  applications  all  distributions  can  be  approximated  with  sufficiant  accuracy 
by  discrete  ones.  Even  for  h  as  small  as  20  the  number  of  terms  in  the  multiple 
sum  computation  gets  easily  out  of  hand  and  the  problem  is  no  longer  practical  to 
solve  by  direct  summation.  This  is  especially  true  because  function  avaluations  are 
computationally  expensive  since  each  term  in  the  multiple  sum  requires  the  solution 
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of  a  linear  program. 


4.  Importance  Sampling 

Monte  Carlo  Methods  are  recommended  to  compute  multiple  integrals  or  mul¬ 
tiple  sums  for  higher  h-dimensional  sample  spaces  (Davis  and  Rabinowitz  (1984), 
Glynn  and  Iglehart  (1989)).  Suppose  Cu  =  C(vu)  are  independent  random  variates 
of  vu,(j  =  l,...,n  with  expectation  z,  where  n  is  the  sample  size.  An  unbiased 
estimator  of  z  with  variance  a\  —  <r2/n,  a2  =  var(C(V))  is 

z  =  (1/n)  £  cr. 

W=1 


Chosing  q*(vu)  =  ^  wou^  ^ea<^  *°  var(*)  —  0>  which  means  one  could 

get  a  perfect  estimate  of  the  multiple  sum  from  only  one  estimation.  Practically 
however,  this  is  useless  since  to  compute  q(vu)  we  have  to  know  z  =  ^wef}  Cup(vu), 
which  we  eventually  want  to  compute.  The  result  however  helps  to  derive  a  heuris¬ 
tic  for  choosing  q.  It  should  be  proportional  to  the  product  C(vu)p(vu)  and  should 
have  a  form  that  can  be  integrated  easily.  Thus  a  function  r(vw)  «  C(vu)  is  sought, 
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which  can  be  integrated  with  less  costs  than  C(vu).  Additive  and  multiplicative 
(in  the  components  of  the  stochastic  vector  v)  approximation  functions  and  combi¬ 
nations  of  these  are  potential  candidates  for  our  approximations.  In  particular,  we 
have  been  getting  good  results  using  C(V)  «  £*=1  C^VJ).  We  compute  q  as 

*  cvwn  ' 

E?=iE 

In  this  case  one  has  to  compute  only  h  1-dimensional  sums  instead  of  1  h-dimensional 
sum.  The  variance  reduction  depends  on  how  well  the  approximation  function  fits 
the  original  cost  function.  If  the  original  cost  function  has  the  property  of  additivity 
(seperability)  the  multiple  sum  can  be  computed  exactly  by  h  1-dimensional  sums. 
If  the  addidive  model  is  a  bad  approximation  of  the  cost  function  the  only  “price” 
that  has  to  be  paid  is  increasing  the  size  of  samples.  If  the  observed  variance  is  too 
high  using  a  starting  sample  size  the  size  of  sample  is  adjusted  higher.  Actually 
we  use  a  variant  of  the  addidive  approximation  function.  By  introducing  C'(r), 
the  costs  of  a  base  case,  we  make  the  model  more  sensitive  to  the  impact  of  the 
stochastic  parameters  v. 

h 

r(V)  =  C(r)  +  £ C(n, .... r,-, ,  Vi, rj+1, . . .  ,n)  -  C(r) 

1=1 

We  denote  this  as  a  marginal  cost  model,  r  can  be  any  arbitrary  chosen  point  of 
the  set  of  values  Vj,  i  =  1, . . . ,  h.  For  example  we  choose  T{  as  that  outcome  of 
which  leads  to  the  lowest  costs,  ceteris  paribus. 

5.  Benders  Decomposition 

We  decompose  the  2-stage  multi-period  stochastic  linear  problem  by  applying 
Benders  (1962)  decomposition.  (Van  Slyke,  Wets  (1969)): 
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The  master  problem: 


T  T 

Zm  =  min^c,  xt  + 

*= 1  1=1 

— Bt-iXt-i  +  At  Xt  —  bf,  t  =  1, ...,T,Bo  —  0 

—G\  r,t  +  ot[0t  >  g[i  t  =  1,  •  •  • ,  T,  /  =  1,  •  •  • ,  L 
xt  >  0 

where  the  latter  constraints,  called  cuts,  are  initially  absent  but  are  inserted  in 
later  iterations.  The  master  problem  is  optimized  to  obtain  an  approximate  optimal 
feasible  solution  x%  =  x\  that  is  used  as  input  to  the  subproblems. 

The  sub-problems  for  uit  in  period  t: 

*"*(*!)  =  min  fty? 

<•(*'):  Dty?'  =d?+F?‘x\,  ut  6  =  1,...,T 

Vt*  >  0,  x\  given , 

where  7r"‘  =  7r"‘(iJ)  are  dual  multipliers  corresponding  to  the  constraints  and 
Zt*  =  is  the  value  of  the  objective  as  functions  of  (x{).  These  are  used  to 

generate  the  next  cut  for  the  master. 

The  cuts:  for  t  =  1, 2, . . . ,  T, 

Gj  =£(*"•  B"'),  ,!  =  £(*”•<!“•),  z,(x[)  =  E{zf),  <•=*,"■(*{). 

Lower  ( LBl )  and  upper  ( UBL )  bounds  to  the  problem: 

T 

LBl  =  zh ,  UBl  =  minttfS1-1,  +  zt(x  j))},  =  oo 

t= l 

x{  is  the  optimal  solution  of  the  master  problem  in  iteration  /,  tt"'(x{)  is  the 
optimal  dual  solution  of  subproblem  ut,  given  x{.  Note  that  if  the  subproblems  are 
infeasible,  a  slightly  different  definition  of  the  cut  is  used,  a  =  0  corresponds  to 
feasibility  cuts  and  a  —  1  to  optimality  cuts.  Solving  the  master  problem  in  iteration 
l  we  obtain  a  trial  solution  x[  which  we  pass  to  the  subproblems.  By  solving  a 


sample  of  sub-problems  u>t,wt  €  Nt,t  —  1  according  to  the  importance 

sampling  scheme  we  compute  estimates  of  the  second  stage  costs  zt  and  estimates 
of  the  gradients  G\  and  the  right  hand  sides  g\  of  the  cuts.  Cuts  represent  in  outer 
linearization  of  the  second  stage  costs  expressed  in  first  stage  decision  variables  and 
6t.  Note  that  there  is  one  cut  for  each  period  t.  The  cuts  are  added  to  the  master 
problem  and  the  master  problem  is  solved  again.  The  objective  function  value  of  the 
master  problem  gives  a  lower  bound  estimate  and  the  total  expected  costs  of  a  trial 
solution  x[,t  =  1, . . .  ,T  gives  an  upper  bound  estimate  to  the  objective  function 
value  of  the  problem.  If  the  lower  and  the  upper  bound  are  sufficiently  close,  which 
is  tested  by  a  Student-t  test,  the  problem  is  considered  to  be  solved.  Lower  and 
upper  bounds  can  be  seen  as  a  sum  of  i.i.d.  random  terms  which  for  sample  sizes 
of  30  or  more  can  be  assumed  normally  distributed  with  known  (derived  from  the 
estimation  process)  variances.  A  95%  confidence  intervall  of  the  optimal  solution 
is  computed.  See  Dantzig  and  Glynn(1990)  and  Infanger(1990)  for  details  of  the 
algorithm. 

6.  Numerical  Results 

This  method  for  solving  large-scale  two-stage  stochastic  linear  programs  with 
recourse  has  been  implemented.  The  code  of  MINOS  (Murtagh  and  Saunders 
(1983))  has  been  adapted  for  this  purpose  as  a  subroutine  for  solving  both  the 
master-problem  and  the  sub-problems.  When  solving  large  numbers  of  sub-problems 
it  is  important  for  the  performance  of  the  algorithm  to  take  advantage  of  good  start¬ 
ing  bases.  Computation  time  can  be  reduced  dramatically  by  solving  first  an  ex¬ 
pected  value  problem  by  replacing  the  stochastic  parameters  by  their  expectations. 
The  expected  value  solution  of  the  resulting  deterministic  problem  is  then  used  as  a 
starting  point  for  the  stochastic  solution.  Additionally  we  keep  cuts  obtained  from 
the  expected  value  problem  to  initially  guide  the  algorithm.  It  can  be  shown  that 
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cuts  obtained  from  the  expected  value  problem  are  valid  for  the  stochastic  problem. 
They  are  “weak”  and  get  repaced  as  the  algorithm  proceeds.  The  code  uses  sparce 
matrix  techniques  and  efficient  data  structures  for  handling  large-scale  problems. 

Computational  results  of  the  large  scale  test  problems  are  represented  in  Ta¬ 
ble  1.  Besides  the  solution  of  the  stochastic  problems,  the  results  from  solving  the 
expected  value  problems  are  also  reported.  We  also  report  on  the  estimated  ex¬ 
pected  costs  if  the  expected  value  solution  is  used  as  the  decision  in  a  stochastic 
environment.  The  objective  function  value  of  the  true  stochastic  solution  has  to  lie 
between  the  minimum  value  of  objective  function  of  the  deterministic  problem  and 
the  expected  costs  of  the  expected  value  solution. 

Expansion  planning  of  multi-area  power  systems 

WRPM  is  a  multi-area  capacity  expansion  planning  problem  for  the  western 
USA  and  Canada.  The  model  is  very  detailed  and  covers  6  regions,  3  demand 
blocks,  2  seasons,  and  several  kinds  of  generation  and  transmission  technologies. 
The  objective  is  to  determine  optimum  discounted  least  cost  levels  of  generation 
and  transmission  facilities  for  each  region  of  the  system  over  time.  The  model  min¬ 
imizes  the  total  discounted  costs  of  supplying  electricity  (investment  and  operating 
costs)  to  meet  the  exogenously  given  demand  subject  to  expansion  and  operating 
constraints.  A  description  of  the  model  can  be  found  in  Dantzig  et.  al.  (1989).  In 
the  stochastic  version  of  the  model  the  availabilities  of  generators  and  transmission 
lines  and  demands  are  subject  to  uncertainty.  There  are  13  stochastic  parame¬ 
ters  per  time  period  (8  stochastic  availabilities  of  generators  and  transmission  lines 
and  5  uncertain  demands)  with  discrete  distributions  with  3  or  4  outcomes.  The 
model  covers  a  time  horizon  of  3  future  periods  of  10  years  each.  Thus  the  total 
number  of  stochastic  parameters  is  39.  The  operating  sub-problems  of  each  period 
are  stochastically  independent.  The  number  of  universe  scenarios  is  larger  than 
5  •  106  per  period.  In  the  deterministic  equivalent  formulation  the  problem  if  it  were 
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possible  to  state  it  would  have  more  than  4.5  billion  constraints. 

The  stochastic  WRPM  is  solved  by  using  a  sample  size  of  100  per  period.  It 
takes  129  iterations  to  obtain  the  expected  value  solution  and  additional  68  it¬ 
erations  to  compute  the  stochastic  solution.  The  objective  function  value  of  the 
stochastic  solution  was  estimated  as  199017.4  with  an  amazingly  small  95%  confi¬ 
dence  intervall  of  0.029%  on  the  left  side  and  0.067%  on  the  right  side.  Thus  the 
optimal  solution  lies  with  95%  confidence  between  198959.3  <  z*  <  199164.1.  The 
expected  costs  of  the  expected  value  solution  (202590.3)  and  the  objective  function 
value  of  the  stochastic  solution  differ  significantly  from  the  expected  costs  of  the 
optimal  stochastic  solution.  The  problem  was  solved  in  687  minutes  on  a  Toshiba 
T5200  laptop  personal  computer.  This  time  includes  time  to  solve  26295  linear 
sub-problems  of  the  size  of  302  rows  and  289  columns  and  197  master  problems. 
Multi-period  Portfolio  Management 

LP42  is  a  portfolio  management  test-problem,  formulated  as  a  network  prob¬ 
lem.  It  is  a  modified  version  of  test-problems  found  in  Mulvey  and  Vladimirou 
(1989).  The  problem  is  to  select  a  portfolio  which  maximizes  expected  returns  in 
future  periods  taking  into  account  the  possibility  of  revising  the  portfolio  in  each 
period.  There  are  also  transaction  costs  and  bounds  on  the  holdings  and  turnovers. 
The  test  problem  covers  a  planning  horizon  of  four  future  periods.  The  returns 
of  the  stocks  in  the  four  future  periods  are  assumed  to  be  independent  stochastic 
parameters,  discretely  distributed  with  3  outcomes  each;  this  formulation  differs 
from  that  of  Mulvey  and  Vladimirou  who  restricted  the  problem  size  by  looking 
at  a  certain  number  of  preselected  scenarios.  Like  in  Mulvey  and  Vladimirou  the 
multi-period  problem  is  viewed  as  a  2-stage  problem,  where  all  future  periods  are  in¬ 
cluded  in  the  second  stage.  With  13  stocks  with  uncertain  returns,  the  problem  has 
52  stochastic  parameters.  The  universe  number  of  scenarios  6  •  1024  is  very  large,  so 
that  the  deterministic  equivalent  formulation  of  the  problem  if  expressed  explicitely 
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would  have  more  than  1.9  •  10 27  rows.  Here,  the  stochastic  parameters  appear  in 
the  B-matrix  as  well  as  in  the  D-matrix.  In  this  case  cuts  from  the  expected  value 
problem  are  not  valid  for  the  stochastic  problem.  The  expected  value  problem  and 
the  stochastic  problem  are  solved  seperately.  A  sample  size  of  600  was  chosen.  The 
solution  (objective  function  value  2.329)  is  obtained  in  4  iterations.  The  95%  confi¬ 
dence  intervaU  is  very  small  given  the  large  number  of  stochastic  parameters,  namely 
0.536%  on  the  left  side  and  0.767%  on  the  right  side.  Thus  with  95%  confidence 
the  objective  function  value  of  the  optimal  solution  lies  within  2.316  <  z*  <  2.347. 
The  expected  costs  of  the  expected  value  solution  is  significantly  different  from  the 
expected  costs  of  the  stochastic  solution. 

7.  Conclusion 

We  have  demonstrated  that  our  approach  using  Benders  decomposition  and 
importance  sampling  is  capable  of  solving  large-scale  problems  of  planning  under 
uncertainty.  Numerical  results  of  large  problems  with  numerous  stochastic  param¬ 
eters  indicate  that  very  accurate  solutions  of  the  problems  can  be  obtained  using 
only  small  sample  sizes. 
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Table  1:  Large-scale  problems:  computational  results 


#  periods 

#  iter  stoch.  (exp.  val.) 
sample  size  (per  period) 
exp.  val.  solution,  obj 
exp.  val.  solution,  exp.  cost 
stochastic  solution,  obj 

est.  conf.  left  % 

conf.  right  % 

solution  time  (min) 


Problem  Size 
Master 


rows 

columns 

nonzeros 

rows 

columns 

nonzeros 


#  stoch.  parameters  (total) 

#  univ.  seen,  (per  period) 


WRPM 

LP42 

3 

1 

197  (129) 

4(6) 

100 

600 

196471.4 

1.611 

202590.3 

2.334 

199017.4 

2.329 

0.0292 

0.536 

0.067 

0.767 

687 

209 

128 

49 

226 

83 

413 

133 

302 

178 

289 

309 

866 

570 

39 

52 

5  •  10® 

6  •  1024 
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